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A theoretical analysis of the thermodynamic response functions of the 2D single-band Hubbard 
model is carried out by means of the composite operator method. It is shown that all the features of 
these quantities can be explained by looking at the dependence of the thermodynamic variables on 
their conjugate ones, that is, for example, the relation between entropy and temperature, chemical 
potential and particle concentration, double occupancy and on-site Coulomb repulsion. In this way, 
the electronic specific heat and the entropy per site are determined in the paramagnetic phase. Also, 
for the electronic specific heat and internal energy we present three different schemes of calculation. 
A comprehensive comparison among them for the interacting case opens the possibility to obtain 
a deep theoretical understanding of the computed quantities whereas, for the non-interacting case, 
only confirms the internal consistency of the thermodynamic approach. It is found that the nu- 
merical data from quantum Monte Carlo techniques for the internal energy and electronic specific 
heat are well reproduced by determining them through the first and second temperature deriva- 
tives of the chemical potential. The anomalous normal state properties in hole-doped cuprate high 
Tc superconductors are also well described. Actually, the results for the linear coefficient of the 
electronic specific heat are in agreement with those obtained by using a pure fermionic theoretical 
scheme. Indeed, the results for the Wilson ratio seem to confirm the same scenario. Finally, we 
obtain several characteristic crossing points for the response functions when plotted versus some 
thermodynamic variables. These peculiar features indicate the existence of more than one energy 
scale competing with thermal excitations and indicate, as already noted by Volhardt, a crossover 
from a non-interacting to a highly correlated behavior. 



I. INTRODUCTION 

It is believed, on both experimental and theoretical grounds, that superconductivity and charge transport in high Tc 
cuprates are mostly confined to the Cu02 planes [1,2]; and hence the attention of many physicists has been dedicated 
to 2D models which contain as an essential feature a competition between the band picture and highly correlated 
many body effects. Of course, some features of the phase diagram, like the existence of a finite Neel temperature, can 
only be explained by adding a coupling between the planes. 

The bonding combination of Cu and O orbitals turns out to be quite deep below the Fermi level, so that no 
dynamical freedom is left to treat d and p orbitals separately [3] (there are some strong experimental evidences, 
mostly based on the study of the Knight shift, that in the Cu02 plane one spin degree of freedom is observed [4]). 
Through the Pauli principle, the energy of the p electron excitation is, for example, largely modified by the change 
of charge and spin states of the neighboring Cu ions. A p electron and charge and spin fiuctuations on neighboring 
Cu ions are simultaneously excited so that electronic excitations are formed on a Cu02 cluster as a whole. Then, the 
resulting complex can be described by a single-band Hubbard model [5] . 

In the simplest form, the Hubbard model, first introduced to describe the correlations of electrons in a narrow d-band 
^ of transition metals, contains a kinetic term which describes the motion of the electrons among the sites of the Bravais 
" " ' lattice and an interaction term between electrons of opposite spin on the same lattice site. By varying the model 
parameters, it is believed that the Hubbard model is capable of describing many properties of strongly correlated 
fermion systems. Among different examples, the Hubbard Hamiltonian is applicable to describe the metal-insulator 
transition in a series of transition metal oxides such as Sri-^LaxTiO^ [6,7] and V2O3 [8-11]. The applicability of the 
model to the superconducting copper-oxides is related to the fact that upon doping most of these compounds exhibit 
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a metal-insulator Mott transition; the superconducting state is near the Ncel state and there are many experimental 
results [12-15] which show a close relation between the antiferromagnetic correlations in the Cu — O planes and the 
occurrence ol the superconducting phase. However, it is important to stress that an appropriate description of a 
bad metal with large energy scale spin fluctuations by means of a purely electrostatic Hamiltonian should preserve 
the symmetry expressed by the Pauli principle that codifies the correct interplay between charge and magnetic 
configurations [16-18]. 

Although considerable attention has been devoted to the Hubbard model and significant progress was achieved in 
understanding ground state properties, particularly at half-filling, static and dynamical spin correlations, the optical 
conductivity and other observables, a clear comprehension of the low- lying excitations is still lacking [2]. The difficulty 
is not to be found only in the absence of any obvious small parameter in the strong coupling regime. More deeply, it is 
due to the difficulty of handling simultaneously itinerant aspects (spatial correlations) and atomic aspects (pronounced 
on- site quantum fluctuations) [19]. 

In recent years we have been developing a method of calculation, denominated Composite Operator Method (COM) 
[17,18,20-32], that has been revealed to be a powerful tool for the description of local and itinerant excitations in 
strongly correlated systems. In previous papers, we considered the Hubbard electronic operators for the determination 
of fundamental excitations. A fully self-consistent calculation of the electronic propagator has been realized by means 
of a constraint with the physical content of the Pauli principle [17, 23-25]. We calculated local quantities, as the 
double occupancy and the magnetic moment [17,23], the energy per site [24], the chemical potential [24], the magnetic 
susceptibility [25, 29], the density of states and the quasi-particle spectra [28, 32]. In all the cases, the results show 
a good agreement with those obtained by numerical simulation. In particular, the results obtained for the magnetic 
properties can reproduce the unusual characteristics observed in high superconducting materials [25, 29, 33]. 
Therefore, the agreement strengthens the idea that a microscopic single-band model contains the essential physical 
features of the new class of materials. 

In this paper we investigate the electronic specific heat and the entropy per site of the 2D Hubbard model for a 
paramagnetic ground state. It will be shown that all the features of these quantities can be understood by looking 
at the dependence of the chemical potential and double occupancy on their conjugate thermodynamic variables, that 
is, the particle concentration and the on- site Coulomb repulsion, respectively. A comprehensive comparison among 
different methods to compute the specific heat will shed a new light on the approximation used. It will emerge that 
in our theoretical scheme, even if dynamical effects in the self-energy are neglected promoting unstable collective 
asymptotic modes to the role of well-defined quasi-particle excitations, extended spin modes can be captured by 
properly combining symmetry requirements and extended operatorial basis. Indeed, the presence of a low temperature 
peak that appears when the low-lying spin states are excited will appear as an important feature shared with the 
quantum Monte Carlo data [34]. In addition, an alternative fermionic scenario for the computation of the electronic 
specific heat together with the results for the Wilson ratio will confirm cuprates as dominated by conventional fermionic 
excitations. An extensive study of the thermodynamic response functions will reveal the existence of critical lines 
which separate different energy scales created by the interplay between charge and spin modes. In other words a 
study of thermodynamics quantities, such as the double occupancy, the entropy, the chemical potential, the specific 
heat, indicates lines in the U — T plane which separate a highly-correlated behavior, dominated by spin and charge 
fluctuations and a non-interacting behavior, dominated by thermal fluctuations. In particular, there emerges a region 
of filling where the entropy reduces by increasing the fllling signalling the set up of an ordered phase. For T ^ 
there is a well-deflned marginal concentration where a quantum phase transition occurs. A detailed comparison with 
the non-interacting case will be also presented throughout the paper. 

The plan of the article is as follows. In the next Section we present the 2D Hubbard model and the electron 
propagator in the COM. In Sec. Ill we review experimental data for some thermodynamic properties. The results 
for the electronic specific heat are presented in Sec IV, where a theoretical understanding of the different ways to 
compute the specific heat is also presented. Section V is devoted to a discussion of double occupancy. In Sec. VI 
the results for the chemical potential versus temperature are discussed. The entropy is analyzed in Sec VII. Some 
concluding remarks are presented at the end. 



II. ELECTRON PROPAGATOR IN THE HUBBARD MODEL 



The Hubbard model is defined by 




(2.1) 
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The notation is the foUowing. The variable i stands for the lattice vector Rj. {c(i), c'l'(i)} are annihilation and 
creation operators of c— electrons at site i, in the spinor notation: 



= (c| c\) (2.2) 



tij denotes the transfer integral and describes hopping between different sites; the U term is the Hubbard interaction 
between two c— electrons at the same site with 

n^{i) = 4(i)c<,(i) (2.3) 

being the charge-density operator per spin a. n{i) is the total charge-density operator, /x is the chemical potential. In 
the nearest neighbor approximation, for a two-dimensional cubic lattice with lattice constant a, we write the hopping 
matrix tij as 



Uj = -Ataij = -4<1 J2 e^'^-^'-'-^j'«(k) (2.4) 

k 

where 



a(k) = ^ [cos(fca;a) + cos(fcj,a)] (2.5) 

The scale of the energy has been fixed in such a way that tu = 0. It should be noted that since the interactions are 
restricted to the same site, the dimensionality of the system comes in only when a specific form for a(k) is taken. In 
other words, the stabilization of eventual cooperative phenomena is uniquely governed by the band dispersion. 

The point of view adopted in the COM is that the Heisenberg operators {c(i),cj(i)} are not good candidates as 
a basis for calculations. Because of strong correlations the c— electrons loose their identity and new fields, whose 
properties are self-consistently determined by the dynamics and by the symmetries of the model, together with the 
boundary conditions, might be more appropriate as a starting point for the physical description of the system. Due to 
the on-site Coulomb interaction, it is known that two sharp features develop in the band structure which correspond 
to the Hubbard subbands and describe interatomic excitations mainly restricted to subsets of the occupancy number. 
Indeed, a first natural choice for composite fields is given by the Hubbard constrained electronic operators 

ai) = [1 - n{iMi) (2.6) 
r,{i) = n{i)c(i) (2.7) 

describing the transitions (n = 0) <;=^ (n = 1) and (n = 1) (n — 2), respectively. The two-point retarded 

thermal Green's function is defined as 

S{i,j) = {R[^(i)^\j)]) , (2.8) 
where ^"(1) is the doublet composite operator 

The bracket (. . .) indicates the thermal average and R is the usual retarded operator. 

In previous papers, we have shown that the determination of the single-particle Green's function (2.8) can be 
realized in a fully self-consistent way once a unique approximation is made [17, 23-25]. This approximation consists 
in neglecting the dynamical part in the self-energy and corresponds to a pole expansion of the spectral intensities. As 
shown in Ref. 23-25, by considering time translational invariance and no magnetic order, the single-particle electronic 
propagator has the following expression 



Scc{i,j) = {RmcHm 



n 



{2n) 



JnB 



Al(k) , A2(k) 



UJ - EiCk) + i-q u) — E2 (k) -|- irj 



(2.10) 



fl and fls being the volume of the unit cell in the direct and reciprocal space, respectively. fF{<^) is the Fermi 
distribution function. £^i(k) and -B2(k) are the energy spectra 
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Ei,2(k) = ii(k)±g(k) 



(2.11) 



Ai(k) and A2(k) are the spectral intensities 

^i(k) = I 



with the definitions 



U{l-n) mi2(k) 
2Q(k) +27n/22Q(k) 



1 + 



U{1 - n) mi2(k) 



2g(k) 27n/22Q(k) 



2 /11/22 



{A + a(k)[p+(l-n)/22]} 



9 \/ - ■ rT72~"''12 
^ V -"11^22 



^ m?,-2[/ii-!^ 

Jll-'22 



TO12 



TOl2(k) = 4t 



A- 



(2.12) 
(2.13) 
(2.14) 

(2.15) 
(2.16) 

(2.17) 

(2.18) 
(2.19) 



where n = (c^c) is the particle density. The parameters A and p are static intersite correlation functions defined as 



(2.20) 



The notation c°'{i) stands to indicate the field c on the first neighbor sites: 



(2.21) 



(2.22) 



We are also using the following notation 



7ii = 1 - n/2 722 = n/2 



(2.23) 



A crucial point in the calculation of the propagator (2.10) is the implementation of the Pauli principle by means 
of constraining equations to determine the internal parameters fi, A and p. Use of this principle leads to the self- 
consistent equations [23, 25] 



n = 1 - Go + C/(l - n)Fo 



(2.24) 



where 



Fn 



pFi = I22F1 - AFo 



Gn 



Br, 



2(27r)2 



Qb 



d^k[a{k)rg{k) 



(2.25) 



2(27r)2 



/ d2fc[a(k)]"/(k)mi2(k) 
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with 



/(k) = 



Ti(k)-T2(k) 
2Q(k) 



5(k) 



[Ti(k)+T2(k)] 



(2.26) 



(2.27) 



T^.(k) 



tanh 



Ei{k) 



2kBT 



) 



(2.28) 



The recovery of the Pauh principle, very often violated by other approximations, assures a dynamics bounded to the 
Hilbert space capable of describing in a correct way the interplay between the charge and the magnetic configurations. 
Furthermore, we have shown [35] that in the two-pole approximation [36] the set of self-consistent equations (2.18) is 
the only one which restores the particlc-holc symmetry and the Pauli principle, which arc intimately connected. 

It is possible to go beyond the two-pole approximation by enlarging the set of asymptotic fields [18, 20-22] or by 
taking into account the dynamical corrections to the self-energy [31,32]. 



The electronic specific heat C{T) of cuprate high-Tg superconductors has been measured. In particular C(T) of 
La^-xSrxCuOi [37, 38] has been studied for 0.03 < x < 0.45 in the range of temperatures between 1.5 and 300 K, 
and of YBa2Cu^Oe+y [39] for 0.16 < y < 0.97 between 1.8 and 300 K. From these experiments the following behavior 

has been observed for the coefficient 7 = C/T of the normal state specific heat: 

a) for fixed temperature, j{x, T) increases with doping; 

al) in the case of La2-xSrxCuOi, 7(3;, T) exhibits a rather sharp maximum at x « 0.25 (near the doping where 
superconductivity disappears), then starts to decrease; the same behavior for 7(0;, T) has been estimated in Ref. 40, 
but with a peak located around x « 0.18, close to the optimal doping; for La2-xBaxCu04 [41] a maximum has been 
observed at a; « 0.22; 

a2) in the case of YBa2CuzO(,+y, j{x, T) increases smoothly to a plateau or two broad maxima, situated at ?/ « 0.6 
and y w 0.9, respectively; 

b) for fixed doping, ^{x,T) as a function of temperature exhibits a broad peak moving to lower temperatures with 
increasing the dopant concentration; 

c) further increasing y, the T-dcpendence weakens and in the region of high doping no increase is observed. For 
Y Ba2Cu30e+y, no substantial increase is observed for y > 0.8. 

As noticed by VoUhardt [42] , there is a peculiar feature of the specific heat observed in a large variety of systems. 
The specific heat curves versus T, when plotted for difi^erent, not too large, values of some thermodynamic variable, 
intersect at one or even two well defined temperatures. In ^He the specific heat C(T, P) curves versus T at different 
pressures P intersect at a well defined temperature [43,44]; in heavy fermions CeAlz [45] and UBeis [46] upon change 
of P, UPt^-xPdx [47] and CePt^Sii-xGex [48] upon change of x, CcCuq-xAux when either P [49] or the magnetic 
field B [50] is varied; in semi-metal, Euq.^Svo.^Ass [51] upon change of B. 

The following properties have been observed for the entropy S [38, 39, 52]: 

a) for a given temperature, S increases with doping; 

al) in the case of La2-xSrxCuOi [38], ^(a;, T) reaches a maximum in the vicinity of a; w 0.25, then decreases; 

a2) in the case of Y Ba2Cu^OQ+y [39, 52], 5* reaches a maximum in the vicinity oi y k, 0.97; 

b) for a given dopant concentration, S exhibits a superlincar dependence on the temperature; 

c) the normal state entropy as a function of T extrapolates to a negative value at T = OJsT; 

d) there is a striking numerical correlation between S/T and axo, where xo is the bulk susceptibility and a is the 
Wilson ratio. 



III. THERMODYNAMICS AS REVEALED BY EXPERIMENTS 



IV. ELECTRONIC SPECIFIC HEAT 



A. GENERAL FORMULAS 



The specific heat C{T) is defined as 



(4.1) 
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where E is the internal energy density, given by the thermal average of the Hamiltonian 



N 



(4.2) 



iV being the number of sites. Calculation of internal energy by means of Eq. (4.2) will generally require the calculation 
of two-particle Green's functions. An alternative way to calculate the internal energy is the following. By introducing 
the Helmholtz free energy per site 



F = E-TS 

where S is the entropy per site, from the thermodynamics we have 

'dF\ fdS 



^J,= 



dn 



dn 



Then, it is straightforward to obtain the following formulas 

rn 

F{T,n)= / iJ.{T,n')dn' 
Jo 



djji 
df 



(4.3) 



(4.4) 



(4.5) 



(4.6) 



E{T,n)= f 
Jo 



KT,n')-T 



df 



dn' 



from which the specific heat turns out to be 



(4.7) 



(4.8) 



In this scheme the thermodynamic quantities are all expressed through the chemical potential, whose determination 
requires the knowledge of the single-particle Green's function. 

Summarizing, we have two distinct ways to calculate the internal energy, based on the use of Eqs. (4.2) and (4.7). 
In principle these equations are equivalent and lead to the same result when an exact solution is available. However, 
the situation drastically changes when approximations are involved and different results can be obtained. Indeed an 
open problem in Condensed Matter Physics is to find a unique consistent scheme of approximation capable of treating 
on an equal footing, both one- and two-particle Green's functions. 

Sometimes in the literature the internal energy for electronic systems is calculated by means of the expression 



E : 



/ + 00 
-OO 



du)N{uj)fF{oj)io 



(4.9) 



where N{u)) is the density of states. However, Eq (4.9) is based on the assumption that the system admits a description 
in terms of fermionic quasi-particles. It is well known that the correlation among the original electrons can generate 
extended bosonic modes, and use of expression (4.9) requires special attention. By general argument it can be shown 
that for interacting systems the correct expression for the internal energy, calculated by means of the density of states, 
is given by 



E: 



/+00 
-oo 



dujN{uj)fF{uj)uJ - {Hi) 



(4.10) 



where {Hi) is the non-quadratic part of the Hamiltonian in the chosen canonical representation. 
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B. NON-INTERACTING CASE 



To discuss the specific heat it is useful at first to consider the non-interacting [i.e. U=0] Hubbard model. In this 
case the thermal retarded Green's function can be exactly calculated and has the expression 

Sco(i,3) = {Rle{i)e<(j)]) 

° L "-^-^^'^ '"-"-'-""-"'li - M^) l _ ^ (4.11) 

where the energy spectrum has the expression 

£;(k) = -/i - 4ia(k) (4.12) 
The chemical potential is determined as a function of n and T by means of the equation 

"=7^12/ d'kME{k)] = l--—^ d'kTik) (4.13) 

where we put 

As we discussed above, we have different ways of calculating the free energy. In the non-interacting case, where 
an exact solution is available, all different procedures must give the same result. Since this point will acquire some 
relevance in the interacting case, we shall examine this in detail. 
In the non-interacting case the density of states is given by 

^^'^^ = (2?F X '^^^^^'^ ' ^^^^^ ^^'^^^ 



By substituting this expression into Eq. (4.10) we have 

Am 



E = 



[ d^ka{k)T{k) (4.16) 
JnB 



where a factor 2 has been included in order to take into account the spin degree of freedom. 
The same result is obtained by taking the thermal average of the Hamiltonian 

E = 8t(c«(i)ct(i)) = -^1 d\aik)T{k) (4.17) 

where use has been made of the expression (4.10) for the single-particle Green's function. 

The proof that use of Eq. (4.7) leads to the same result requires some work. By taking the derivative with respect 
to T of Eq. (4.12) we have 



^^I _ 1 

df ~ T 



(4.18) 



where 



K = A / i"^'^^" (4-19) 



By considering that the derivative with respect to n of Eq. (4.12) gives 

V, Ja^d'kaikf-^ n d 



substitution of (4.17) into (4.8) leads to 



d^ka(k)T{k) (4.20) 
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FIG. 2. 7(r, n) of the non-interacting 2D Hubbard model is plotted as a function of the filling for various values of the 
temperature. 



C. INTERACTING CASE 



In the interacting case the expressions (4.2) and (4.7) for the internal energy give different results. By recalling 
(2.10), Eq. (4.2) gives 

Eh =Ef- UD (4.23) 

where D is the double occupancy. We use Eh to indicate that the internal energy has been calculated by the average 
of the Hamiltonian. Ep is the part of the internal energy calculated on the assumption of only fermionic elementary 
excitations 

/ + 0O 
duiN{u})fF{u})u (4.24) 
-OO 

= (^X d'H{^-Ti{k)][E,ik)+n]A,{k) + [l-T2{k)][E2{k)+n]A2ik)} 

(4.23) 

Alternatively, wc have seen that the internal energy and the specific heat can be calculated by means of Eqs. 
(4.8) and (4.9). This procedure requires a knowledge of the first and second temperature derivatives of the chemical 
potential. Let us define 

By taking the first derivative with respect to T of the self-consistent equations (2.18), we obtain 

At = ^G(^) - ) + l-^B[^^ (4.26) 

Z Z ZI11I22 



where 

9T™ 9T™ gj^m \^-^'J 

Explicit calculation of the derivatives defined in Eq. (4.26) shows that the equations (4.25) provide a set of linear 
algebraic equations for the three parameters /Ui, Ai, pi, as functions of the parameters fx, A, p. In the same way, by 
taking the second derivative with respect to T of Eqs. (2.18) we obtain 

G^^^ = U{l-n)F^'^ 

A2 = l^Gf ) - ) + ) (4.28) 



P2F1 + 2piF[^^ +pF[^^ = /22Ff ) - A2F0 - 2AiFo(^) - AFq^'^ 

These equations provide a set of linear algebraic equations for the three parameters /i2, A2, P2 as functions of the 
parameters /U, A, p, m, Ai, pi. Once the self-consistent calculation of the three parameters /x. A, p has been performed 
by means of the set (2.18), then the calculation of the first and second derivatives of the chemical potential reduces 
to the solution of simple linear equations. 

In the interacting case, because of the approximation used, the different procedures to calculate the internal energy 
give different results. At first we shall compare our theoretical results with the data obtained by numerical analysis. 
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principles. 

The paper was devoted to the study of the electronic specific heat and entropy per site in the paramagnetic phase. 
We analyzed these quantities by looking at the dependence of the thermodynamic variables on their conjugate ones, 
that is, for example, the relation between entropy and temperature, chemical potential and particle concentration, 
double occupancy and on-site Coulomb repulsion. Once the self-consistent equations for the single-particle propagator 
have been solved, we have determined the temperature derivatives of the internal parameters by means of exact linear 
systems of algebraic equations. The determination of the first and second temperature derivatives of the chemical 
potential has been revealed crucial in determining the thermodynamic response functions under investigation. For 
the electronic specific heat and internal energy we have presented three different schemes of calculation. All of them 
allowed the possibility to obtain a deep theoretical understanding of how and to which extent collective excitations 
can be retained in the description of thermal response functions. We have obtained a good agreement with the 
data by quantum Monte Carlo techniques for the electronic specific heat and the internal energy [34]. Further on, 
although Monte Carlo data shared common features with the results from the calculations through the T-derivative 
of the chemical potential, the experimental data for cuprates, as revealed by the Wilson ratio and linear coeSicient 
of the electronic specific heat, have shown that in such systems the dominant excitations are those of conventional 
non-interacting fermions [38] . 

We obtained several characteristic crossing points for the response functions when reported as functions of some 
thermodynamic variables. These peculiar features, already evidenced by Vollhardt [42], marked turning points where 
diff'erent response functions evolve from a non-interacting behaviour 

(i) the entropy is an increasing function of U; 

(ii) the entropy is an increasing function of n; 

(iii) the double occupancy is a decreasing function of T; 

(iv) the T-derivative of the chemical potential is a decreasing function of n; 

(v) the linear coefficient of the specific heat is an increasing function of n; 
to an unconventional dependence on the conjugate variables 

(vi) the entropy is a decreasing function of U; 

(vii) the entropy is a decreasing function of n; 

(viii) the double occupancy is an increasing function of T; 

(ix) the T-derivative of the chemical potential is an increasing function of n; 

(x) the linear coefficient of the specific heat is a decreasing function of n. 

Before closing wc would like to mention that the region of filling, where (vi)-(x) hold, coincides with that where 
incommensurate magnetism and superconductivity are experimentally observed in LSCO cuprates family. 
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